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Abstract. The ab initio theory of electronic excitations in atomically thin [quasi- 
two-dimensional (Q2D)] crystals presents extra challenges in comparison to both the 
bulk and purely 2D cases. We argue that the conventionally used energy-loss function 
—Im (where e, q, and w are the dielectric function, the momentum, and the 

energy transfer, respectively) is not, generally speaking, the suitable quantity for the 
interpretation of the electron-energy loss spectroscopy (EELS) in the Q2D case, and 
we construct different functions pertinent to the EELS experiments on Q2D crystals. 
Secondly, we emphasize the importance and develop a convenient procedure of the 
elimination of the spurious inter-layer interaction inherent to the use of the 3D super¬ 
cell method for the calculation of excitations in Q2D crystals. Thirdly, we resolve 
the existing controversy in the interpretation of the so-called ir and n + a excitations 
in monolayer graphene by demonstrating that both dispersive collective excitations 
(plasmons) and non-dispersive single-particle (inter-band) transitions fall in the same 
energy ranges, where they strongly influence each other. 


1. Introduction 

Quasi two-dimensional (Q2D) crystals (the systems periodic and infinite in two 
dimensions while of atomic-size thickness) are presently attracting much of attention 
because of their potential in future nano-technologies. In particular, the electronic 
excitations in Q2D crystals, such as electron-hole pairs generation, plasmons, excitons, 
etc., are of great importance, since they determine the response to external actions, 
and, ultimately, shape the materials’ functionality (see, e.g., [I] for a recent review). At 
the same time, the theoretical approaches as well as the computational methods to deal 
with the excitations in Q2D systems become much more involved compared to both the 
three-dimensional (3D) and purely 2D cases, and are still far from finally established. 

Indeed, in both purely 3D and purely 2D cases (the latter referring to model systems 
of zero thickness), one usually introduces the dielectric function e(q, cu), where q is the 
wave-vector of the corresponding dimensionality, and oj is the frequency. The definition 
of e reads [j] 

= £{q,U})(j) tot {q,uj), (1) 

where 4> ext and (j) tot are the externally applied and the total (external plus induced) 
scalar potentials of the electric field, respectively. The definition ([TJ) is not, however, 
straightforwardly transferable to the Q2D case, which can be easily appreciated by 
considering that (f) tot depends on the z (perpendicular to the layer) coordinate, even 
when (p ext is uniform in z direction. Therefore, we need to redefine the dielectric 
function (energy-loss function, conductivity, etc.) of a Q2D crystal in accordance to 
the system’s structure. Importantly, the new definitions must be consistent with the 
quantities measured experimentally, the violation of which requirement having caused 
much of confusion in recent literature. 

f In this paper we are concerned with the longitudinal fields. 
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Secondly, the widest used at this time method of dealing with Q2D systems 
numerically is the super-cell geometry approach. This is to artificially replicate the 
system periodically in z direction, choosing the distance d between the layers large 
enough to prevent the inter-layer interactions. The resulting fictitious system is 3D- 
periodic, and well developed 3D solid-state methods as well as the existing 3D software 
can be used to efficiently handle it. In the ground-state calculations, the interpretation 
of the super-cell results is straightforward: Provided that the inter-layer distance is such 
that the corresponding wave-functions do not overlap, properties of an isolated single 
layer coincide with those of a single super-cell. 

However, if we are concerned with the dynamic response to electric held, the 
situation changes. In the case of a Q2D crystal, the electric held decays into vacuum as 
~ e~ q \ z \. Clearly, whatever large is the separation d between the layers, at sufficiently 
small q the inter-layer interaction persists, and results of the super-cell calculation cannot 
be literally transferred to a single layer. A simple method of taking d very large does 
not work, because small q-s are often of special interest, and the computational load 
grows rapidly with the increase of d. Therefore, we need to establish relations between 
the response functions of the hctitious 3D array system and those of the single layer 
of interest. The paramount requirement to these relations (test for sanity) is that they 
should produce the same response functions of a single layer from the differing 3D 
response functions calculated with different separations d. 

In this paper, in the framework of the time-dependent density-functional theory 
(TDDFT), we construct a scheme satisfying the above requirements. As its practical 
application, we calculate and perform a detailed analysis of the excitation spectrum of 
a single-layer graphene in the 1-30 eV energy range. Resolving the recent controversy 
in the interpretation of the tt and it + a peaks as plasmons la si m E] or single-particle 
inter-band transitions [6], we show that, while the problem is far more complicated 
than thought before, there do exist prominent tt and 7r + cr plasmons in graphene, 
although they overlap with the inter-band transitions. The isolation of the both kinds 
of excitations becomes possible by following the wave-vector dependence of spectra 
(spatial dispersion). We also show that much of the previous confusion was due to 
using unsuitable theoretical quantities for comparison with experiment, as well as to 
the uncritical use of the results of the super-cell calculations in application to Q2D 
systems. 

This paper is organized as follows. In Sec. [2]we re-examine and revise the definitions 
of the response-functions (dielectric function, energy-loss function, conductivity) 
pertinent to measurements on Q2D crystals. Sec. [3] is devoted to the methods of 
obtaining the response-functions of Q2D crystals from the 3D super-cell calculations. 
In Sec. [4] we present results of calculations performed for graphene as an important 
example of a Q2D crystal, discuss the differences arising from the use of our theory, 
and compare with experiment. Section [5] contains conclusions, and in the Appendix 
we give the derivation of some formulas presented in Sec. [3] We use the atomic units 
(e 2 = h = m e = 1) unless explicitly specified otherwise. 
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Quasi-2D crystals are neither periodic three-dimensionally nor strictly two-dimensional 
within the mathematical plane: They rather have a finite atomic or nano-scale thickness. 
Accordingly, one has to exercise caution in the choice of quantities (dielectric function, 
conductivity, density-response function, etc.) which characterize the dielectric response, 
giving them rigorous, suitable for Q2D systems, meaning. The safe way to do this 
is to start from the general 3D case without assuming the periodicity (translational 
invariance) in any direction. The density-response function y(r, r', oj) is defined in the 
real-space representation as 

p{r,w) = J x( r , r ', uj)(j) ext (r', u>)dr', (2) 

where (p ext is the scalar potential of the externally applied electric field and p is the 
charge-density induced in the system in response to (p ext . The inverse generalized 
dielectric function e _1 (r, r', to) is defined as 

= J £■ 1 (^,r / ,n;)0 e2:^ (r , ,a;)dr , , (3) 

where 0 is the total (external plus induced) potential, and it is expressed through x as 

£ _1 (r, r', u) = 8{v - r') + [ (4) 

J | r — r | 

By Fourier-transform, Eqs. (j2])-(j4j) can be equivalently written in the reciprocal-space 

representation, of which we will need only 

An t 

£_1 (q, q r , u ) = 5(q- q') + — x(q,q , ,w)- (5) 

qz 

For Q2D crystals, we will use the mixed reciprocal- and real-space representation, 
in the xy plane and z direction, respectively. Then the definition (j2j) becomes 

OO 

PGn(*,q||,w) =Y^ j (6) 

f-1 / ’’ 

°ll -oo 

where Gy and Gy are the 2D reciprocal lattice vectors, and qy is the 2D wave-vector in 
the first Brillouin zone. We will also need the expression for the induced potential 

OO 

2t r f 

<(*-. qi |,w) = |G J p G ,(z',<l h u)e-' a ^«‘-‘''dz'. (7) 

—OO 

We distinguish between different kinds of experimental setups, where it will be 
shown that the measurable quantities correspond to different definitions of the energy- 
loss function. 

2.1. Electron energy-loss spectroscopy 

The principal experimental method of probing qy- and ^-dependent response of Q2D 
systems is the Electron Energy-Loss Spectroscopy (EELS). We will consider separately 
the reflection and transmission EELS. 
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2.1.1. Reflection EELS. In reflection EELS, the quantity, which characterizes the 

Its definition reads: Let 


inelastic electron scattering, is the so-called (/-function pfi 
the external potential bctjj] 

4> ext (z, qy, <n) = e qz . 


( 8 ) 


Then asymptotically at z —> +oo 

fto d (z,<l\\,u) = g{q\\,uj)e- qz . ( 9 ) 

With the use of Eqs. (EJ) and ([7]), the (/-function can be written as [TI] [T2, [7] 

OO 

g(q\\,uj) = — J e q ^ z+z,) x 00 (z,z',q h uj)dzdz'. (10) 

—OO 

When qn and u> are substituted with the incident electron’s loss of the parallel 
momentum and energy, respectively, the differential cross-section of the inelastic 
electron scattering in the reflection geometry is determined by the energy-loss function 
-Im#(q||,w), i.e., 

d 2 cr in (p' p) . . . . 

—— -L(q ( , W ) = -Imj(q,, W ), (11) 

where dQ is the element of the solid-angle of the detection of the scattered electron, 
q = p — p', lo = E — E', E = p 2 / 2, E' = p' 2 / 2, p and p' are the momenta of the 
incident and scattered electron, respectively. Equation (fTUj) expresses g through the 
density-response function y. 


2.1.2. Transmission EELS. For the transmission EELS of an arbitrary (possibly, 
non-periodic) target, the differential cross-section of the inelastic electron scattering 
is proportional to the energy-loss function (see, e.g., Pi) 


d 2 <J in { p' 1- p) 


U<1| :CJi --tIim f (q,q,u 


( 12 ) 


dLldu ' 

We emphasize, that in the right-hand side of Eq. (TT2T) of Eq. ([5]) enters with q and 
q' both substituted with the same momentum transfer q. 

We consider the case of normal incidence. Then 


( P - 9.) 2 , 9 _ 

+ --E-U, 


Qz — P 

and, since p is large, 


p- 


2uj 


qf. 


Qz = 


2u + q 2 ^ 

2 p 


(13) 

(14) 

(15) 


§ In this paper, we restrict ourselves to the dipole scattering regime [7]. Impact-scattering regime is 
more involved, requiring the solution of the inelastic scattering problem within the distorted-waves 
method ®mm, which is outside the scope of this paper. 

|| Everywhere below, only the z dependence of quantities is written down explicitly. The implied ry 
and t dependence is always e^ q ll' r ll _wt h 
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Then, by Eqs. (IT!?|) and (151) , and restricting ourselves to the momentum transfer qy 
within the first Brillouin zone, we have 

OO 

L(q,u) = -^ImJ e iqz{z '~ z) x 00 {z,z',q\\,u)dzdz'. (16) 

— OO 

Equation (fTTfl) is the valid energy-loss function to be used in the theory of the 
transmission EELS. We point out that it replaces —^ Im^d__ |14j . the latter not being 
defined for Q2D crystal. 


2.1.3. In-plane electronic transport. Let us calculate the conductivity and the energy 
dissipation in a Q2D crystal under the action of the uniform in z direction external 
potential 


fi ex \z, q||,(n) = 1. 

Then, by Eq. ([6]), 

OO 

Po(z,q\\,uj) = J Xoo{z,z\q l{ ,uj)dz', 

— OO 

and for the 2D particle-density averaged in the xy plane 

OO OO 

p2D{q\\,u)= j Po(z,q\\,u)dz = J Xoo{z,z',qu,uj)dzdz' 

—OO —OO 

Using the continuity equation 

Q|| ’ j 2 o(q||,w) = ujp 2D (q\\,u), 

where j 2 d is the 2D current-density, we can write 


OO 


janlqi.w) = J Xoo(z,z',q l ,u)dzdz', 


(17) 


(18) 


(19) 


( 20 ) 


( 21 ) 


and since, by Eq. ffTTjl . E ext = — iqy, 

j 2 Z?(q||,w) = a™£(q\\,u)E ext (q\\,u), (22) 

where a is the 2D conductivity with respect to the external field, which is given by 

OO 

= "T / Xoo(z,z',q h uj)dzdz'. (23) 

q w L 

For the energy dissipation Q per unit time per unit cell area under the action of 
the external potential (fI71) we can writc^ 


27 T / UJ OO 

<3(q|h") = 2 ^ J dt J dz J *|j(r,«)-E(r,«), 


(24) 


0 —oo A 


Because Q is the quadratic rather than linear property, we must temporarily return to the (r, t) 
representation. 
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where j(r, £) and E(r ,t) are well defined 3D current-density and electric field, 
respectively, and A is the 2D unit cell of the Q2D crystal. Since carriers do not commit 
work on themselves, Eq. ()24l) can be rewritten as 

2n/ui oo 

<3( q| b^) = ^ J dt J dz J dry j(r,t) ■ E ext (r,t) (25) 

0 — oo A 


and simplified to 


Q(qib^) = o Re 


E ex (q||,w)* • / i 0 (z,q\\,u)dz 


= -^Re 
2 q 2 


E'(q||,w)*-q|| / p 0 (z,q\\,u)dz 


(26) 


Using Eqs. (J6J) and (fT71) . we can rewrite Eq. (1261) as 


<5(q|b w ) = yRecx^qy,^ = Im / Xoo(z,z , ,q\\,u)dzdz'. 


CO 


(27) 


The important conclusion of Sections 12.1.21 - 12.1.31 is that the energy losses 
in the reflection EELS, transmission EELS, and the energy dissipation in the 
in-plane transport are governed by the quantities J e q ^ z+z lxoo( z , z', qy, uo)dzdz', 
f e t9z ^ z '~ z ^xoo(z, z', qy, co)dzdz', and f Xoo(z,z',q\\,oo)dzdz', respectively, which, for a 
Q2D crystal, are clearly different quantities. We note that in the limiting case of the 
purely 2D (zero-thickness) crystal, x( z , U, qy, oo) contains 8(z)8(z'), and, consequently, 
all these three quantities coincide. 


2.1.4- Dielectric function of Q2D crystal. It has already been noted in Sec. [I] that 
the dielectric function cannot be rigorously introduced for a Q2D crystal in a simple 
multiplicative way via Eq. fill) . This is, in the first place, due to the uncertainty at 
which z f> tot (z, qy,w) must be used. It is, however, convenient to keep the notion of 
the dielectric function on the intuitive level, resorting to an analogy with the purely 2D 
system. I 11 the latter case 


Xoo(z,z',q h u) = 8(z)8(z')x oo(q||,w), ( 28 ) 

where x^P G , (qy, co) is the 2D density-response function and the relation holds 

1 0 77- 

—7-T = 1 + —Xoo(qib^), (29) 

C 2 D(q||,w) q\\ 

where c 2 d( q, w) is the dielectric function of the 2D system. By Eqs. (1231) and 


ICO 


= — Xoo(qii^) 

q \\ 


a 2D (qib 


(30) 
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and, using Eq. 


= 1 + 


27r 9|l 


6 2D { q||,w) tuj 


-a. 


2D 


(qih")- 


(31) 


For a Q2D crystal, we define the dielectric function using Eq. ([3T]) , where in the 
right-hand side we substitute erfof Eq. (123|) . which is a well defined quantity for the 
Q2D crystal. Therefore 


£Q2d(<1\\, w) 


= 1 + 


2vr 




(32) 


By comparing Eqs. (ITUll . (HU . and (1321) . we see that the inelastic scattering of 

electrons in EELS of Q2D system is not defined by the loss-function —Im- } --. It 

will be shown that this difference is of a quantitative importance in the interpretation 
of the measurements on Q2D crystals. 

Below we will see that it is convenient to have our results in the full reciprocal-space 
representation. Using the Fourier expansion in the [—§,!] interval, where at this point 
we consider d an arbitrary sufficiently large distance, we can write 

SG z z 


e kz = 


^2 b G z (k)t 


(33) 


G z 


where 


h G z {k) = 


2 sinh 


(. k-iG z )d 
2 


d{k — iG z ) 
Then Eq. (ITU]) can be written as 

47t d 


L{ q,w) = 


Eq. m as 


and Eq. 


T 


2nd 


Im E b G z z) b G' z (iQz)X0G z ,0G , z (q||, w) 


G Z G' 


f/(q||,w) = - b GM\) b G' z (q\\)XOG z ,OG'M\\,u), 


as 


911 


. , . 

— 19 -Xoo.oofqihw) 


(34) 


(35) 


(36) 


(37) 


e Q 2D(q||,w) q\\ 

Having established expressions for the observables via the density-response function 
of a single layer x G||G ' (u E, qy, a;), we turn to the methods of the calculation of the 
latter. 
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= 11111 




Figure 1. Schematics of 2D material under an external field, a) 2D single-layer 
geometry and b) 3D super-cell geometry. 


3. Excitation of Q2D system from the super-cell geometry calculation 


By far the widest used at this time way of the calculation of the electronic response of 
atomically thin Q2D crystals is the super-cell geometry approach. This is, instead of 
considering a single layer [schematized in Fig. [lj a)], to artificially periodically replicate 
it in z-direct ion [Fig. |T] b)], and study the resulting 3D periodic system in place of 
the original single-layer one. It is, however, known puna mi that the 3D dielectric 
function obtained in the super-cell geometry has nothing to do with that of a single¬ 
layer system of interest, unless a special procedure of extracting the 2D response from 
the 3D super-cell geometry calculation (i.e., the elimination of the spurious inter-layer 
interaction) is applied. In Appendix A we give a detailed derivation of the expression for 
the density-response function y of the single-layer system from y - the density-response 
function of the array system. It reads in the matrix formE] 


x(qib^) = x(qii^) [i + ^(qii)x(qii,^)] 

where the matrix Cqg' is given by 


-i 


(38) 


CGG'(qil) - F Gz G' z (|G|| + q|||)5 G ||G' 
4vr (p 2 - G Z G' Z ) 


G Z G‘ 


(p) = 


pd{p 2 + G 2 z )(p 2 + G'l) 


cos 


(G~ + G' z )d 


(l-e~ pd ). 


(39) 

(40) 


The density-response function of the array system y G|| G' (q||,^) can be routinely 
obtained with the 3D solid-state periodic codes such, e.g., as Elk, abinit, Wicn2k, etc.. 
While y is different for different layers’ separations d, y obtained through Eq. (j38l) does 
not depend on d, which is crucial for its being a true characteristic of a stand-alone 


+ In Eq. (1381) . the full 3D reciprocal space representation for both y and y is used, which is natural 
for the former, but may look artificial for the latter. We, however, note that an arbitrary function of 
z and z' can be expanded in the Fourier series in z,z' € [— |, |], and then, if necessary, the real-space 
representation can be retrieved by the inverse Fourier transform. 
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2D systenfl We note that Eqs. (I5SD - (TUI are valid not only in RPA [IB], but also in 
TDDFT with (semi-)local f xc (such as, e.g., ALDA), when different layers interact by 
the Coulomb potential only. 

We conclude this section by noting that the ignoring the fundamental difference 
between the single-layer response of a Q2D crystal and that of the array of those layers 
has led to the misinterpretation of the 7r and n + a peaks as single-particle inter-band 
transitions rather than plasmons [6]. 

4. Results of calculations and discussion 



0 5 10 15 20 25 

u (eV) 

Figure 2. Dielectric function of the pristine monolayer graphene obtained with 
Eq. (1371) and the corresponding energy-loss function. Momentum q|| is varied along 
the TM direction. 

We have conducted the super-cell geometry calculation for the monolayer pristine 
graphene followed by the application of Eq. / fdffi) for the extraction of the single-layer 
response (elimination of the spurious inter-layer interaction) from the 3D calculation. In 
the DFT super-cell calculation, we used the full-potential linear augmented plane-wave 

* Strictly speaking, this is XG||G' (z, z', qy, ui) which is independent on d. Its Fourier transform in the 
interval z,z' £ |] does formally depend on d, which is of no physical consequence. 
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0 5 10 15 20 25 

u (eV) 

Figure 3. The EELS-related energy-loss function — Im g(q||, w) of Eq. (l36l) (upper 
panel) compared with the energy-loss function —Im fQ2£) ( q °f Eq. (l37l) (lower panel) 
for pristine monolayer graphene. 

(FP-LAPW) code Elk [T9]. The 2 -axis period of the super-cell was taken 20 a.u. The 
fc-point grid of 512 x 512 x 1, 30 empty bands, and the damping parameter of 0.002 a.u. 
were used in both the ground-state and the linear-response calculations. The former was 
conducted within the local-density approximation for the exchange-correlation potential, 
while the latter was the random-phase approximation (RPA) one (i.e., the exchange- 
correlation kernel f xc was set to zero). The dielectric matrix of the 3D system was of 
the size 55 x 55, which was inverted to obtain £ 3 D(q||,^; d) with the local-field effects 
included both in the perpendicular and in-plane directions. 

Results for the Q2D dielectric and energy-loss functions obtained through Eq. (1371) 
are presented in Fig. [2] Features at the energy-loss function around 5 and 15 eV are 
usually referred to as 7r and tt + a peaks, respectively, [2, [3j 0] [6, 5] (lower panel of 
Fig, [2D- To determine whether they are collective excitations (plasmons) [2, 3J 0] [5j or 
single-particle (inter-band) transitions [6], we scrutinize the dielectric function. 

A clear criterion for an excitation to be classified as plasmon is the requirement 
of the real-part of the dielectric function to cross zero at the corresponding energy. 
Otherwise, if the peak at the energy-loss function is due to a peak at the imaginary 
part of the dielectric function, the excitation is a single-particle transition. In our case, 
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Monolayer, n Bilayer, n 



Monolayer, 7r + a 


Bilayer, n + cr 




10 12 14 16 18 20 22 24 12 14 16 18 20 22 24 

u (eV) u (eV) 


Figure 4. The reflection EELS energy-loss function — Im< 7 (q||, u) of Eq. (f36l> (solid 
lines) and the energy-loss function —Im eQ2D ( q|| UJ ) of Eq. (1371) (dashed lines) for 
monolayer (left) and bilayer (right) graphene in the energy-range of 7 r (upper) and 
7 r + <7 (lower) plasmons. Symbols are experimental reflection EELS of Ref. [18] . The 
theoretical spectra have been roughly normalized to the experimental ones. 


starting from greater wave-vectors (from q fa 0.05 and 0.1 a.u., for n and n + a peaks, 
respectively) the real part of the dielectric f un ction does cross zero, and, thereby, from 
those q-s up the peaks are plasmons unambiguously. Are they also plasmons at smaller 
q- s, or they turn into single-particle transitions, although keeping full similarity to their 
higher-g plasmon counterparts ? To answer this question, we need to establish whether, 
although Re C 2 d does not cross zero at smaller wave-vectors, the features at the energy- 
loss function are due to Re 62.0 approaching zero, or they are due to peaks at Ime 2 D (cf., 

m- 

As can be seen in Fig. [2l middle panel, Im £ 2 d has prominent non-dispersive 
(standing in place with the variation of gy) feature at oj ~ 4 eV. A prominent dispersive 
peak in the energy-loss function (bottom panel) follows the position of not only zero (at 
greater values of q») but also the minimum (at smaller q\\-s) of Re e 2 d hi the range of 4 
- 6 eV, and must, therefore, in both cases be classified as n plasmon. 

The situation is similar, though less transparent, in the case of the 7 r + a peak. 
Here, the broad feature at L(qy,a;), which spreads from u> fa 11 eV upward, is due 
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Figure 5. Dispersion of 7 r and ir + cr plasmons. Momentum q|| is changed along the 
TM direction. 

to the overlap of the single-particle and collective excitations. Up to 00 ~ 14 eV, the 
excitation is purely due to the single-particle transitions and is non-dispersive. Then, 
from 14 eV upward, a dispersive feature due to zero or small values of Re E 2 D begins. 
Because of Re £ 2 d having a low slope and remaining small in absolute value in a wide 
u range, this feature becomes very broad. It, however, is clearly plasmon as well. 

Interestingly, in two points, at lo ~ 8 and 14 eV, Re £20 becomes nearly unity for all 
values of q\\, so that all the curves in Fig. [2j upper panel, intersect at these two points. 
The physical reason for this is still to be understood. 

I 11 Fig. [3] the reflection EELS-related energy-loss function — Img(q||,a;), calculated 
by Eq. (l36lh is compared with the loss-function —Im £Q2D ( q|| of Eq. (I37lh the latter 
relevant to the energy dissipation in the in-plane transport . While the 7r features (~ 
5 eV) are similar for both loss-functions, the 7 r + a features are considerably different, 
the EELS loss-function having a richer structure in this energy range. This can be 
understood by noting that, according to Eq. (1371) . 1/e is determined by the head element 
of the density-response matrix Xoooo(q||j w), while, by Eq. (l36lh the ^-function also 
includes contributions from XoG z OG'„( f l||; w) with nonzero G z and/or G' z . The role of 
the latter increases with the increase of c 0 , transitions into the higher bands becoming 
allowed, which leads to the differences in the spectra. 

In Fig. [4] we compare our theory with reflection EELS experiment of Ref. [18] on 
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monolayer and bilayer graphene. In the energy range of 7r plasmon (two upper graphs 
in Fig. H]), ^-function and the energy-loss (L) function are practically identical and both 
compare to experiment rather well. On the contrary, in the range of 7T + a plasmons, the 
g- and L-functions acquire differences for both monolayer and, particularly, for bilayer 
graphene (lower graphs). For bilayer graphene 7r + cr peak at (/-function is narrower than 
at the L- function, the former being closer to experiment (lower right). Remarkably, for 
monolayer graphene both g- and L-functions find themselves in shear disagreement with 
experiment (lower left). 

Obviously, the main source of the disagreement between the theory and experiment 
can be the presence of the SiC (0001) substrate in the latter and the absence thereof 
in the former. The substrate may cause additional screening leading to the shift of 
plasmon peaks and change the dispersion law [20] j2L, ;22j [23]. Moreover, the interaction 
of graphene with a substrate can cause the deformation of the graphene sheet, leading 
to the confinement of the plasmon oscillations between the ripples [24] , 

Let us briefly discuss the modifications to the theory necessary to include Q2D 
crystals supported on substrates. For the ab initio theoretical treatment of surfaces, 
either pure or adsorbates covered, the super-cell method is widely applied as well. 
In variance with the stand-alone Q2D systems, for surfaces, rather thick slabs of the 
material alternated with the similar thickness vacuum slabs are used. Equations (1381) - 
(1401) are readily applicable in this case too, producing as an output the density-response 
function of a single slab. This slab having two surfaces is not, however, the semi-infinite 
system of interest. Neither the array of such slabs, from which we have started, is. Which 
system, the array of slabs or a single slab, is better to model a Q2D crystal on top of a 
semi-infinite substrate must be decided in each particular case. An alternative method 
may be the inclusion of the substrate by means of the effective background dielectric 
function after having calculated the response of a stand-alone Q2D system. Although it 
may be useful in some situations, this would necessarily be a phenomenological approach, 
probably inconsistent with the rigorousness of the calculation for the Q2D crystal alone. 
Finally, the ultimate theory will be that asymptotically matching Q2D crystal on the 
surface with the 3D periodic bulk of the substrate, which is a very challenging task. 

The neglect of the many-body exchange-correlation effects [25, [26], both static and 
dynamic, may further contribute to the discrepancy. Both the accurate inclusion of the 
substrate in the theoretical setup and the inclusion of the many-body effects by means 
of using nonlocal exchange-correlation kernel of TDDFT remain the major challenges 
of the present-day theory of Q2D materials. 

Figure [5] shows the dispersion of tt and 7r + cr plasmons with the qy vector varied 
in the TM direction calculated in the EELS setup. The 7r + cr plasmon has a complex 
structure and is split into distinct sub-excitations. 

As noted in Ref. [6], in a 2D system there can be no plasmon at gy = 0. In the 
ordinary 2D electron gas [27] and in doped graphene in the case of 2D plasmon due 
to electrons in the n* band [28j, this is ensured by the plasmon disappearance via its 
energy tending to zero as y/gj[ when gy —> 0. However, as can be seen from Fig. [2] lower 
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panel, for n and 7r + cr plasmons in graphene the same realizes in a different way: The 
energy positions of these plasmons converge to finite values (4 and 14 eV for n and n + cr 
plasmons, respectively) as q\\ —* 0, while the amplitudes of the corresponding peaks go 
to zero. In this regard, as in many others, graphene is also special compared to the 
ordinary 2D electron gas. The zero limit of the plasmons’ amplitude can be understood 
from simple arguments. Indeed, quite generally, in a 2D crystal, the relation between 
the dielectric function e 2D and the conductivity cr 2 D is 

^(q,.^) = l + 2 ^" <T2g(q "- Cj) . (41) 

CO 

Since in the single-layer graphene, a 2 D{<i\\ = 0,(u) is finite [and known [29j to tend to 
e 2 /4:h as u —* 0], £ 2 £>(q|| = 0,o;) is unity identically, leading to the zero energy-loss 
function. Our calculated conductivity of the monolayer pristine graphene at q» = 0 
versus the frequency is shown in Fig. 0 where the static limit of 1/4 is shown by the 
dotted line. 



to (eV) 


Figure 6. Frequency-dependent conductivity of pristine monolayer graphene at 
q\\ = 0. When u —> 0, the conductivity must converge to its analytic two-bands 
model value of 1/4, which it does, except for too small u> (at the very tip of the Dirac 
cone), where the accuracy of the calculation is limited by the /c-point grid. 
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We have identified and dealt with the extra challenges the ab initio theory of electronic 
excitations in quasi-2D crystals presents in comparison to both the bulk and purely 2D 
cases. In particular, we have re-examined the problem of the correspondence between 
the theoretically calculated quantities and the observables in the measurements on quasi- 
2D crystals and found that the energy-loss function —Im conventionally used for 

the interpretation of the EELS data, is not, generally speaking, the right quantity to 
be compared with this kind of the experiment. Instead, the EELS-related energy-loss 
functions, specific for both transmission and reflection experimental setups, must be 
used, which has been shown to be both qualitatively and quantitatively different in the 
case of quasi-2D systems. In the limit of the zero crystal thickness, our theory reduces 
to the conventionally used one. 

We have addressed the problem of the classification of the tt and 7r + a peaks in 
monolayer graphene. By the use of the accurate procedure of extracting the dielectric 
function of a 2D crystal from the 3D super-cell geometry calculation, we have shown 
conclusively that there exist prominent 7r and 7r + a collective exitations (plasmons) in 
graphene, although they are accompanied by interband transitions in the close energy 
ranges. Plasmons and single-particle interband transitions can be distinguished from 
each other provided the wave-vector-resolved analysis is performed. 

We have also demonstrated the importance of the correct interpretation of the 
super-cell geometry calculation results, the latter taken simplistically, can lead to 
erroneous qualitative conclusions for materials of reduced dimensionality. We expect 
that the present theory of the electronic excitations in Q2D crystals, by correctly 
accounting for the atomic or nano- thickness of the systems studied, will further boost 
the theoretical and experimental work, improving the methods of comparison between 
the two. 
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Appendix A. Extraction of the response of a single layer from that of the 
array-of-layers system 

Let us consider a periodic array of identical layers, as schematized in Fig. [U and let 
us construct the density-response function y G||G / (z, z> qy, w) of the array system. We 
assume that the separation d between the layers is large enough so that densities, both 
the ground-state and the dynamically induced, do not overlap between layers. Then for 
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d/2 


n G 




G' 


II ~d/2 


where 


<t> G f {z,q.\\,u) = #f(*) 
2vr 


d/2 


|G + q||| 

1 N 1 m —— 


E' 


n G (z', qy,cu)e _|G+q||ll ^ _2,_md| (iE. (A.2) 


—d/2 


Equation (1A.ll) uses the fact that the potential external to the rn — 0 layer is that of 
Eq. (1A.2I) . i.e., it is the proper external potential plus the potentials induced by all the 
layers with m^O. The prime at the sum means that the m — 0 term is excluded from 
the summation. Performing the summation explicitly, we hav^l 

</4 f (z,q h u) =<l>G t {z,m\,u) 

d/2 

4-7T f 

n G {z' , q||, oj) cosh[|G + qy | [z — z')\dz'. (A.3) 


G + qi||(el G+q iil d -l) 


-d/2 

Expanding into the Fourier series 

cosh [p(z - z')\ = d g z G' z (p) ( 


,iG z z—iG',z' 


GzG' 


where 


n , 4(p 2 - G Z G' Z ) 

D G z G'{P) = -rror COS 


d?(p 2 + G 2 z )(p 2 + G' z 
and introducing the matrices 


(G z + G' z )d 


sinli 


2 I P d 


(A.4) 


(A.5) 


F Gz G'(p ) = 


And 


p ( eP d _ ^ D GzG’ z (p) 

4vr(p 2 - G Z G' Z ) 


cos 


(G z + G' z )d 


P d (p 2 + G 2 )(p 2 + G' z 2 ) 

C'G||G*,G'Gi(q||) = F Gz G' z ( |G|| + q|| |)<5g,|G' , 
we arrive at the fully 3D reciprocal-space representation 

^(qipw) = (j> e G (q.\\,uj) + Y C G G '(<l\\)n G '(q.\\,u), 

G' 

«G(q||,w) = 


(1 - e~ pd ) 


(A.6) 

(A.7) 

(A.8) 

(A.9) 


jj Equation (IA.3I) corrects an error in Eq. (A3) of Ref. [T7] . However, all the formulas derived in 
Ref. |17j under the assumption of qa <C 1, where a is the effective width of a layer, remain valid. 
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Using Eqs. (1A.8j> and (lA.9p . the definition of \ as the density-response function of 
the array system 

n G (q||,a;) = ^XGG'(q||, w)0^(q||, u), (A. 10 ) 

G' 

and in view of the arbitrariness of (p ext , we can write in the matrix form 


X(q||>w) = x(q||,w) + x(q|hw)Cx(q||,w). 


(A.ll) 


Finally, inverting Eq. (1A.11IL we arrive at Eq. (J38]). 

It may be useful to consider approximations to the above exact scheme. First, since 
G||d 1 for G|| 7 ^ 0, this is practically exact to keep in Eq. (1A.2I) the Gy = 0 term 
only. Further, if qa 1, where a is the effective width of the layer, Eq. (I38]i can be 
shown to yield a convenient relation between the corresponding 3D and 2D dielectric 
functions P2! 


l 

£Q 2 D(q||,w) 



1 


1 


1 


. 

s3D(q||> w ; 

d) 

911“ 


+ 



(A.12) 


where £ 3 £>(q||, d) is the dielectric function of the fictitious periodic 3D system 
comprised of the layers separated by the distance d. 

In the q\\d 1 limit we have from Eq. (IA.12D 


= i + ‘Ml 


£2D(qn,w) 2 L e 3r>(qn,w;d) 

In the long-wave limit q\\d <C 1 

_±_ y = i + ^ [ i- £sD(qil , w;d) ], 

which, since the second term is small, can be written as 

£2D(q||.w) = 1- [1 -£3d(<J||,W;<0] . 


(A.13) 


(A.14) 


(A.15) 


[1] Politano A and Chiarello G 2014 Nanoscale 6(19) 10927 10940 

[2] Gass M H, Bangert U, Bleloch A L, Wang P, Nair R R and Geim A K 2008 Nature Nanotechnology 

3 676-681 

[3] Eberlein T, Bangert U, Nair R R, Jones R, Gass M, Bleloch A L, Novoselov K S, Geim A and 

Briddon P R 2008 Phys. Rev. B 77(23) 233406 

[4] Liu Y, Willis R F, Emtsev K V and Seyller T 2008 Phys. Rev. B 78(20) 201403 

[5] Liou S C, Shie C S, Chen C H, Breitwieser R, Pai W W, Guo G Y and Chu M W 2015 Phys. Rev. 

B 91(4) 045418 

[6] Nelson F J, Idrobo J C, Fite J D, Miskovic Z L, Pennycook S J, Pantelides S T, Lee J U and 

Diebold A C 2014 Nano Letters 14 3827-3831 

[7] Liebsch A 1997 Electronic excitations at metal surfaces (New-York: Plenum) 

[8] Nazarov V U 1995 Surf. Sci. 331/333 1157 

[9] Nazarov V U 1999 Phys. Rev. B 59 9866 

[10] Nazarov V and Nishigaki S 2001 Surface Science 482 - 485 640 647 

[11] Persson B N J and Zaremba E 1985 Phys. Rev. B 31 1863 

























Electronic excitations in quasi-2D crystals 


19 


[12] Tsuei K D, Plummer E, Liebsch A, Pchlkc E, Kempa K and Bakshi P 1991 Surface Science 247 

302 - 326 

[13] Sturm K 1982 Adv. Phys. 31 1 

[14] Wachsmuth P, Hambach R, Kinyanjui M K, Guzzo M, Benner G and Kaiser U 2013 Phys. Rev. B 

88(7) 075433 

[15] Rozzi C A, Varsano D, Marini A, Gross E K U and Rubio A 2006 Phys. Rev. B 73(20) 205119 

[16] Despoja V, Novko D, Dekanic K, Sunjic M and Marusic L 2013 Phys. Rev. B 87(7) 075447 

[17] Nazarov V U, Alharbi F, Fisher T S and Kais S 2014 Phys. Rev. B 89(19) 195423 

[18] Lu J, Loh K P, Huang H, Chen W and Wee ATS 2009 Phys. Rev. B 80(11) 113410 

[19] http://elk.sourceforge.net 

[20] Politano A, Marino A R, Formoso V, Farias D, Miranda R and Chiarello G 2011 Phys. Rev. B 

84(3) 033401 

[21] Politano A, Marino A, Formoso V, Faras D, Miranda R and Chiarello G 2012 Plasmonics 7 369-376 

[22] Generalov A and Dedkov Y 2012 Carbon 50 183 - 191 

[23] Cupolillo A, Politano A, Ligato N, Perez D C, Chiarello G and Caputi L 2015 Surface Science 

634 76 - 80 

[24] Politano A, Campi D, Formoso V and Chiarello G 2013 Phys. Chem. Chem. Phys. 15(27) 11356- 

11361 

[25] Botti S, Sottile F, Vast N, Olevano V, Reining L, Weissker H C, Rubio A, Onida G, Del Sole R 

and Godby R W 2004 Phys. Rev. B 69 155112 

[26] Nazarov V U and Vignale G 2011 Phys. Rev. Lett. 107(21) 216402 

[27] Stern F 1967 Phys. Rev. Lett. 18 546-548 

[28] Das Sarma S and Hwang E H 2009 Phys. Rev. Lett. 102(20) 206412 

[29] Wunsch B, Stauber T, Sols F and Guinea F 2006 New Journal of Physics 8 318 



